
### Set working directory to current folder ###
# if sourcing in base R
tryCatch({
  setwd(getSrcDirectory(function(){})[1])
}, error=function(e){})

# if running in RStudio
tryCatch({
  setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
}, error=function(e){})


### Initialize relevant packages ###
library(haven)
library(dplyr)
library(tidyr)
library(stringr)
library(fixest)
library(ggplot2)
library(kableExtra)

### Initialization ###

## Read in data

load("TaxCompliance_Data_Lab.RData")
data_conjoint <- read_dta("Conjoint.dta")
data_lab_conjoint <- read_dta("Lab_Conjoint.dta")

### MAIN RESULTS, LAB-IN-THE-FIELD ###

## Figure 3 ## 
ggplot(data_lab, aes(x = round_order, y = tax)) +
  geom_point(stat = "summary", fun = "mean", size = 2.5,
             color = "black", alpha = .9) +
  geom_line(stat = "summary", fun = "mean", linetype = "dashed") +
  geom_vline(xintercept = 7) +
  geom_line(data = data_lab %>%
              group_by(ST_treat_group, ST_post) %>%
              mutate(tax_mean = mean(tax)) %>%
              group_by(ST_treat_group) %>%
              mutate(tax_mean_pre = ifelse(round_order <= 7,
                                           max(tax_mean * (1 - ST_post)),
                                           NA),
                     tax_mean_post = ifelse(round_order >= 7,
                                            max(tax_mean * (ST_post)),
                                            NA)),
            aes(x = round_order,
                y = tax_mean_pre,
                group = ST_treat_group),
            size = 1.25, alpha = .75) +
  geom_line(data = data_lab %>%
              group_by(ST_treat_group, ST_post) %>%
              mutate(tax_mean = mean(tax)) %>%
              group_by(ST_treat_group) %>%
              mutate(tax_mean_pre = ifelse(round_order <= 7,
                                           max(tax_mean * (1 - ST_post)),
                                           NA),
                     tax_mean_post = ifelse(round_order >= 7,
                                            max(tax_mean * (ST_post)),
                                            NA)),
            aes(x = round_order,
                y = tax_mean_post,
                group = ST_treat_group),
            size = 1.25, alpha = .75) +
  facet_wrap(~ ST_treat_group) +
  ggtitle("Tax Compliance Rate by Experiment Round") +
  xlab("Round") +
  ylab("Tax Compliance Rate") +
  scale_x_continuous(breaks = c(1:12),
                     labels = c(1:12)) +
  scale_y_continuous(labels = function(x) sub("0\\.", "\\.", x)) +
  coord_cartesian(ylim = c(.5,1)) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 11.5),
        axis.text.y = element_text(angle =90,
                                   vjust = .5, hjust = .5,
                                   size = 11.5),
        plot.title = element_text(face = "bold", size = 12),
        axis.title = element_text(size = 12),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank())
ggsave("fig_3.png",
       width = 6.75, height = 5, units = "in",
       dpi = 300)

## Table 2 ##
list(
  # immediate effect, no FE
  feols(tax ~ treat, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(round_order %in% c(6,8)),
        cluster = ~ respondent),
  # immediate effect, respondent FE
  feols(tax ~ treat | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(round_order %in% c(6,8)),
        cluster = ~ respondent),
  # average effect, respondent FE
  feols(tax ~ treat | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned),
        cluster = ~ respondent),
  # average effect, respondent FE, common time trend
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend
  feols(tax ~ treat | respondent[round_order], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent)
) %>% etable(digits = "r3",
             digits.stats = "r3",
             tex = TRUE,
             title = "Impact of different types of government officials on
             tax compliance",
             fitstat = c("n", "r2"),
             keep = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             order = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             dict = c(`tax`="Tax Compliance",
                      `treatNon-Responsive`="Non-Responsive",
                      `treatResponsive`="Responsive",
                      `treatNon-Sanctioning`="Non-Sanctioning",
                      `treatSanctioning`="Sanctioning",
                      `respondent`="Respondent"),
             drop.section = "slopes",
             file = "table_2.tex",
             replace = TRUE,
             label = "tab:result_itt",
             float = TRUE,
             headers = list("Immediate Effect"= 2, "Average Effect" = 3),
             extralines=list("-\\emph{Time trends}"=c("", "", "", "", ""),
                             "-Common"=c("", "", "", "Yes", ""),
                             "-Individual"=c("", "", "", "", "Yes")),
             placement = "htbp!")

### APPENDICES ###

## Table A1 ##
summary_table_conjoint <- data_conjoint %>%
  # clean up data and do calculations
  filter(task == 1) %>%
  dplyr::select(age,
                male,
                edu,
                ccp,
                local_hukou,
                income_cat,
                apartment_owned,
                married) %>%
  mutate_at(vars(age:married), as.integer) %>%
  summarise_at(vars(-group_cols()),
               list(mean = ~mean(., na.rm = TRUE),
                    sd = ~sd(., na.rm = TRUE),
                    n = ~sum(!is.na(.)))) %>%
  pivot_longer(age_mean:married_n,
               names_to = c("variable", "measure"),
               names_pattern = "(.+)_(.+)") %>%
  pivot_wider(names_from = measure,
              values_from = value) %>%
  # transform into table
  mutate(mean = formatC(mean, digits = 3, format = "f"),
         sd = paste("(", formatC(sd, digits = 3, format = "f"), ")",
                    sep = "")) %>%
  mutate(variable = case_when(variable == "age" ~ "Age (Numeric)",
                              variable == "male" ~ "Gender (Binary, Male = 1)",
                              variable == "edu" ~ "Education (Ordinal, 1-6)",
                              variable == "ccp" ~ "CCP (Binary)",
                              variable == "local_hukou" ~ "Hukou (Binary)",
                              variable == "income_cat" ~ "Income (Ordinal, 1-23)",
                              variable == "apartment_owned" ~ "\\# of Apartments Owned (Ordinal, 1-5)",
                              variable == "married" ~ "Married (Binary)",
                              TRUE ~ variable)) %>%
  mutate(variable = paste(str_replace(variable, "\\(", "{\\\\itshape \\("), "}"))


sink("table_a1.tex")
kable(summary_table_conjoint,
      format = "latex",
      align = "lccc",
      booktabs = TRUE,
      escape = FALSE,
      row.names = FALSE,
      caption = "Summary Statistics for Conjoint Experiment Respondents",
      label = "summary_stats_conjoint",
      col.names = c("", "Mean", "Std. Dev.", "Number of Respondents"),
      linesep = "") %>%
  kable_styling(latex_options = c("hold_position", "scale_down"),
                position = "c")
sink()

## Table A2
## Balance Table for SSI 
feols(c(age, male, edu, ccp, local_hukou, income_cat, apartment_owned, married) ~ 
        exemption + penalty_3 + citizen_input + prov_supervise,
      data = data_conjoint %>% filter(source == 0),
      cluster = ~respondentindex) %>% 
  etable(digits = "r2",
         digits.stats = "r2",
         tex = TRUE,
         title = "Balance Test for SSI Conjoint Experiment Data",
         fitstat = c("n", "r2"),
         order = c("Exemption", "Penalty", "Citizen Input", "Top-down Supervision"),
         dict = c(`age`="Age",
                  `male`="Gender",
                  `edu`="Education",
                  `ccp`="CCP",
                  `local_hukou`="Hukou",
                  `income_cat` = "Income",
                  `apartment_owned` = "\\# Apts. Owned",
                  `married` = "Married",
                  `exemption`="Exemption",
                  `penalty_3`="Penalty",
                  `citizen_input`="Citizen Input",
                  `prov_supervise`="Top-down Supervision",
                  `respondentindex`="Respondent"),
         file = "table_a2.tex",
         replace = TRUE,
         label = "tab:balance_ssi",
         float = TRUE,
         placement = "htbp!",
         adjustbox = 1)

## Table A3
## Balance Table for Qualtrics 
feols(c(age, male, edu, ccp, local_hukou, income_cat, apartment_owned, married) ~ 
        exemption + penalty_3 + citizen_input + prov_supervise,
      data = data_conjoint %>% filter(source == 1),
      cluster = ~respondentindex) %>% 
  etable(digits = "r2",
         digits.stats = "r2",
         tex = TRUE,
         title = "Balance Test for Qualtrics Conjoint Experiment Data",
         fitstat = c("n", "r2"),
         order = c("Exemption", "Penalty", "Citizen Input", "Top-down Supervision"),
         dict = c(`age`="Age",
                  `male`="Gender",
                  `edu`="Education",
                  `ccp`="CCP",
                  `local_hukou`="Hukou",
                  `income_cat` = "Income",
                  `apartment_owned` = "\\# Apts. Owned",
                  `married` = "Married",
                  `exemption`="Exemption",
                  `penalty_3`="Penalty",
                  `citizen_input`="Citizen Input",
                  `prov_supervise`="Top-down Supervision",
                  `respondentindex`="Respondent"),
         file = "table_a3.tex",
         replace = TRUE,
         label = "tab:balance_qualtrics",
         float = TRUE,
         placement = "htbp!",
         adjustbox = 1)



## Table C4 ##
## Summary Table
summary_table_lab <- data_lab %>%
  # clean up data and do calculations
  filter(round_order == 1) %>%
  dplyr::select(ST_treat_group,
                age_group,
                male,
                edu,
                ccp,
                local_hukou,
                income_cat,
                apartment_owned,
                married) %>%
  mutate_at(vars(age_group:married), as.integer) %>%
  group_by(ST_treat_group) %>%
  summarise_at(vars(-group_cols()),
               list(n = ~n(),
                    mean = ~mean(., na.rm = TRUE),
                    sd = ~sd(., na.rm = TRUE))) %>%
  mutate(n = age_group_n) %>%
  dplyr::select(-ends_with("_n")) %>%
  pivot_longer(age_group_mean:married_sd,
               names_to = c("variable", "measure"),
               names_pattern = "(.+)_(.+)") %>%
  pivot_wider(names_from = measure,
              values_from = value) %>%
  mutate(se = sd/sqrt(n)) %>%
  dplyr::select(ST_treat_group, n, variable, mean, se) %>%
  # transform into table
  mutate(mean = formatC(mean, digits = 3, format = "f"),
         se = paste("(", formatC(se, digits = 3, format = "f"), ")",
                    sep = "")) %>%
  pivot_longer(mean:se,
               names_to = "measure",
               values_to = "value") %>%
  dplyr::select(ST_treat_group, n, variable, value, measure) %>%
  pivot_wider(names_from = c("ST_treat_group", "n"),
              values_from = "value") %>%
  dplyr::select(-measure) %>%
  mutate(variable = case_when(variable == "age_group" ~ "Age (Ordinal, 1-5)",
                              variable == "male" ~ "Gender (Binary, Male = 1)",
                              variable == "edu" ~ "Education (Ordinal, 1-6)",
                              variable == "ccp" ~ "CCP (Binary)",
                              variable == "local_hukou" ~ "Hukou (Binary)",
                              variable == "income_cat" ~ "Income (Ordinal, 1-23)",
                              variable == "apartment_owned" ~ "\\# of Apartments Owned (Ordinal, 1-5)",
                              variable == "married" ~ "Married (Binary)",
                              TRUE ~ variable)) %>%
  mutate(variable = ifelse(row_number() %% 2 == 0,
                           paste("\\hspace{10pt}", "{\\itshape", str_match(variable, "\\(.*\\)"), "}", "\\vspace{3pt}"),
                           str_remove_all(variable, "\\(.*\\)"))) %>%
  rbind(.,
        c("", rep(NULL,4)),
        c("Obs.",
          str_extract_all(colnames(.),
                          pattern = "\\d+",
                          simplify = TRUE)[-1] %>%
            as.numeric))

sink("table_c4.tex")
kable(summary_table_lab,
      format = "latex",
      align = "lcccc",
      booktabs = TRUE,
      escape = FALSE,
      row.names = FALSE,
      caption = "Summary Statistics by Treatment Conditions",
      label = "summary_stats_lab",
      col.names = NULL,
      linesep = "") %>%
  add_header_above(c("",
                     "Non-Responsive",
                     "Responsive",
                     "Non-Sanctioning",
                     "Sanctioning"),
                   align = "c") %>%
  kable_styling(latex_options = c("hold_position", "scale_down"),
                position = "c")
sink()

## Table C6 ##
## Drop Round 1 and 12

list(
  # average effect, respondent FE
  feols(tax ~ treat | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned) %>%
          filter(round_order > 1, round_order < 12),
        cluster = ~ respondent),
  # average effect, respondent FE, common time trend
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(round_order > 1, round_order < 12),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend
  feols(tax ~ treat | respondent[round_order], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(round_order > 1, round_order < 12),
        cluster = ~ respondent)
) %>% etable(digits = "r3",
             digits.stats = "r3",
             tex = TRUE,
             title = "Impact of different types of government officials on
             tax compliance, excluding first and last rounds",
             fitstat = c("n", "r2"),
             keep = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             order = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             dict = c(`tax`="Tax Compliance",
                      `treatNon-Responsive`="Non-Responsive",
                      `treatResponsive`="Responsive",
                      `treatNon-Sanctioning`="Non-Sanctioning",
                      `treatSanctioning`="Sanctioning",
                      `respondent`="Respondent"),
             drop.section = "slopes",
             file = "table_c5.tex",
             replace = TRUE,
             label = "tab:result_drop_first_last",
             float = TRUE,
             headers = list("Average Effect" = 3),
             #group = ,
             extralines=list("-\\emph{Time trends}"=c("", "", ""),
                             "-Common"=c("", "Yes", ""),
                             "-Individual"=c("", "", "Yes")),
             #highlight = .("gray!80,se"="Sanctioning"),
             placement = "htbp!",
             notes = "Immediate effect estimates are not estimated because
             they do not make use of data from the first and last rounds.")

## Table C6 ##
## Drop Round 7

list(
  # average effect, respondent FE
  feols(tax ~ treat | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned) %>%
          filter(round_order != 7),
        cluster = ~ respondent),
  # average effect, respondent FE, common time trend
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(round_order != 7),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend
  feols(tax ~ treat | respondent[round_order], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(round_order != 7),
        cluster = ~ respondent)
) %>% etable(digits = "r3",
             digits.stats = "r3",
             tex = TRUE,
             title = "Impact of different types of government officials on
             tax compliance, excluding seventh round",
             fitstat = c("n", "r2"),
             keep = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             order = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             dict = c(`tax`="Tax Compliance",
                      `treatNon-Responsive`="Non-Responsive",
                      `treatResponsive`="Responsive",
                      `treatNon-Sanctioning`="Non-Sanctioning",
                      `treatSanctioning`="Sanctioning",
                      `respondent`="Respondent"),
             drop.section = "slopes",
             file = "table_c6.tex",
             replace = TRUE,
             label = "tab:result_drop_7",
             float = TRUE,
             headers = list("Average Effect" = 3),
             #group = ,
             extralines=list("-\\emph{Time trends}"=c("", "", ""),
                             "-Common"=c("", "Yes", ""),
                             "-Individual"=c("", "", "Yes")),
             #highlight = .("gray!80,se"="Sanctioning"),
             placement = "htbp!",
             notes = "Immediate effect estimates are not estimated because
             they do not make use of data from the 7th round.")

## Table C7 ##
## Session FEs

list(
  # immediate effect, session FE
  feols(tax ~ treat | data_session, data = data_lab %>%
          mutate(treat = ST_treat_assigned) %>%
          filter(round_order %in% c(6,8)),
        cluster = ~ data_session),
  # average effect, session FE
  feols(tax ~ treat | data_session, data = data_lab %>%
          mutate(treat = ST_treat_assigned),
        cluster = ~ data_session),
  # average effect, session FE, common time trend
  feols(tax ~ treat + round_order | data_session, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ data_session),
  # average effect, session FE, session time trend
  feols(tax ~ treat | data_session[round_order], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ data_session)
) %>% etable(digits = "r3",
             digits.stats = "r3",
             tex = TRUE,
             title = "Impact of different types of government officials on
             tax compliance, session fixed-effects",
             fitstat = c("n", "r2"),
             keep = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             order = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             dict = c(`tax`="Tax Compliance",
                      `treatNon-Responsive`="Non-Responsive",
                      `treatResponsive`="Responsive",
                      `treatNon-Sanctioning`="Non-Sanctioning",
                      `treatSanctioning`="Sanctioning",
                      `respondent`="Respondent",
                      `data_session`="Session"),
             drop.section = "slopes",
             file = "table_c7.tex",
             replace = TRUE,
             label = "tab:result_session_FE",
             float = TRUE,
             headers = list("Immediate Effect"= 1, "Average Effect" = 3),
             #group = ,
             extralines=list("-\\emph{Time trends}"=c("", "", "", ""),
                             "-Common"=c("", "", "Yes", ""),
                             "-Session"=c("", "", "", "Yes")),
             #highlight = .("gray!80,se"="Sanctioning"),
             placement = "htbp!")

## Table C8 ##
## Alternative time trends

list(
  # average effect, respondent FE, common time trend, linear
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, common time trend, quadratic
  feols(tax ~ treat + round_order + round_order^2 | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, common time trend, cubic
  feols(tax ~ treat + round_order + round_order^2 + round_order^3 | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, common time trend, cubic
  feols(tax ~ treat + round_order + round_order^2 + round_order^3 + round_order^4 | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend, linear
  feols(tax ~ treat | respondent[round_order], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend, quadratic
  feols(tax ~ treat | respondent[round_order] + respondent[round_order^2], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend, cubic
  feols(tax ~ treat | respondent[round_order] + respondent[round_order^2] + respondent[round_order^3], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # average effect, respondent FE, individual time trend, cubic
  feols(tax ~ treat | respondent[round_order] + respondent[round_order^2] + respondent[round_order^3] + respondent[round_order^4], data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent)
) %>% etable(digits = "r3",
             digits.stats = "r3",
             tex = TRUE,
             title = "Impact of different types of government officials on
             tax compliance, alternative time trends",
             fitstat = c("n", "r2"),
             keep = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             order = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             dict = c(`tax`="Tax Compliance",
                      `treatNon-Responsive`="Non-Responsive",
                      `treatResponsive`="Responsive",
                      `treatNon-Sanctioning`="Non-Sanctioning",
                      `treatSanctioning`="Sanctioning",
                      `respondent`="Respondent"),
             drop.section = "slopes",
             file = "table_c8.tex",
             replace = TRUE,
             label = "tab:result_trends",
             float = TRUE,
             headers = list("Average Effect" = 8),
             #group = ,
             extralines=list("-\\emph{Time trends}"=c("", "", "", "", "", "", "", ""),
                             "-Common"=c("Linear", "Quadratic", "Cubic", "Quartic", "", "", "", ""),
                             "-Respondent"=c("", "", "", "", "Linear", "Quadratic", "Cubic", "Quartic")),
             #highlight = .("gray!80,se"="Sanctioning"),
             placement = "htbp!")

## Figure C5 ##
## Event History Plots

event_base <- function(main){
  plot(1,
       xlim = c(1,12),
       ylim = c(-.5, .5),
       main = main,
       xlab = "Round",
       ylab = "Estimated Treatment Effect",
       axes = FALSE,
       pch = 16)
  axis(side = 1, at = seq(1, 12, 1), cex.axis=.8)
  axis(side = 2, at = seq(-.5,.5,.25), labels = c("-.5", "-.25", "0", ".25", ".5"))
  box()
  abline(h = seq(-.5, .5, .25), lty = 2, col = "grey")
  abline(h = 0, lty = 1, col = "black")
  abline(v = 7, lty = 5, col = "black")
}

png("fig_c5_a.png",
    width = 4.25, height = 3, units = "in", res = 300, pointsize = 11.5)
event_base("Non-Responsive")
feols(tax ~ i(round_order, treat, 7) | respondent,
      data = data_lab %>%
        mutate(treat = (ST_treat_group == "Non-Responsive")) %>%
        filter(treat)) %>%
  iplot(ci.width = 0,
        pt.join = TRUE,
        add = TRUE)
dev.off()

png("fig_c5_b.png",
    width = 4.25, height = 3, units = "in", res = 300, pointsize = 11.5)
event_base("Responsive")
feols(tax ~ i(round_order, treat, 7) | respondent,
      data = data_lab %>%
        mutate(treat = ST_treat_group == "Responsive") %>%
        filter(ST_treat_assigned %in% c("Control", "Responsive"))) %>%
  iplot(ci.width = 0,
        pt.join = TRUE,
        add = TRUE)
dev.off()

png("fig_c5_c.png",
    width = 4.25, height = 3, units = "in", res = 300, pointsize = 11.5)
event_base("Non-Sanctioning")
feols(tax ~ i(round_order, treat, 7) | respondent,
      data = data_lab %>%
        mutate(treat = ST_treat_group == "Non-Sanctioning") %>%
        filter(ST_treat_assigned %in% c("Control", "Non-Sanctioning"))) %>%
  iplot(ci.width = 0,
        pt.join = TRUE,
        add = TRUE)
dev.off()

png("fig_c5_d.png",
    width = 4.25, height = 3, units = "in", res = 300, pointsize = 11.5)
event_base("Sanctioning")
feols(tax ~ i(round_order, treat, 7) | respondent,
      data = data_lab %>%
        mutate(treat = ST_treat_group == "Sanctioning") %>%
        filter(treat)) %>%
  iplot(ci.width = 0,
        pt.join = TRUE,
        add = TRUE)
dev.off()

## Table C9 ##
## Heterogeneity By # of Homes Owned and Property Value

list(
  # Model 4 in main table
  # average effect, respondent FE, common time trend
  # entire sample of homeowners
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag),
        cluster = ~ respondent),
  # one home only
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(apartment_owned == 1),
        cluster = ~ respondent),
  # second home and above
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(apartment_owned > 1),
        cluster = ~ respondent),
  # property value in bottom 75%
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(property_value < 8000000),
        cluster = ~ respondent),
  # property value in top 25%
  feols(tax ~ treat + round_order | respondent, data = data_lab %>%
          mutate(treat = ST_treat_assigned_lag) %>%
          filter(property_value > 8000000),
        cluster = ~ respondent)
) %>% etable(digits = "r3",
             digits.stats = "r3",
             tex = TRUE,
             title = "Impact of different types of government officials on
             tax compliance by property ownership and value",
             fitstat = c("n", "r2"),
             keep = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             order = c("Non-Responsive", "Responsive", "Non-Sanctioning", "Sanctioning"),
             dict = c(`tax`="Tax Compliance",
                      `treatNon-Responsive`="Non-Responsive",
                      `treatResponsive`="Responsive",
                      `treatNon-Sanctioning`="Non-Sanctioning",
                      `treatSanctioning`="Sanctioning",
                      `respondent`="Respondent"),
             drop.section = "slopes",
             #file = "table_c9.tex",
             replace = TRUE,
             label = "tab:result_property",
             float = TRUE,
             headers = list(":_:" = list("Full Sample" = 1,
                                         "\\# Apts. Owned" = 2,
                                         "Property Value" = 2),
                            " " = list(" " = 1,
                                       "1" = 1,
                                       "2+" = 1,
                                       "Bottom 75%" = 1,
                                       "Top 25%" = 1)),
             #group = ,
             extralines=list("-\\emph{Time trends}"=c("", "", "", "", ""),
                             "-Common"=c("Yes", "Yes", "Yes", "Yes", "Yes"),
                             "-Individual"=c("", "", "", "", "")),
             #highlight = .("gray!80,se"="Sanctioning"),
             placement = "htbp!")
